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In the high-temperature phase of QCD, the heavy quark momentum diffusion constant 
determines, via a fluctuation-dissipation relation, how fast a heavy quark kinetically equi- 
librates. This transport coefficient can be extracted from thermal correlators via a Kubo 
formula. We present a lattice calculation of the relevant Euclidean correlators in the gluon 
plasma, based on a recent formulation of the problem in heavy-quark effective field theory 

(HQET). We find a « 20% enhancement of the Euclidean correlator at maximal time sep- 
O ■ 

aration as the temperature is lowered from 6T C to 2T C , pointing to stronger interactions at 
■ lower temperatures. At the same time, the correlator becomes flatter from 6T C down to 

' 2T C , indicating a relative shift of the spectral weight to lower frequencies. A recent next- 

to-leading order perturbative calculation of the correlator agrees with the time dependence 
of the lattice data at the few-percent level. We estimate how much additional contribution 
from the ui < T region of the perturbative spectral function would be required to bring it in 
agreement with the lattice data at 3.1T C . 

<D '. 

^ PACS numbers: 12.38. Gc, 12.38. Mh, 25.75.-q 



I. INTRODUCTION 
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CNJ ' Imagine a heavy quark moving through the plasma of light quarks and gluons with more than its 

fair share of thermal kinetic energy f^ksT. Interactions with the medium will gradually cause it to 
slow down, at a rate determined by the drag coefficient rj. It also experiences stochastic interactions 
which balance the drag force in such a way that the heavy quark's kinetic energy tends to |/cbT. 
In this article we study from first principles the strength of these stochastic interactions, which is 
characterized by a parameter called the momentum diffusion coefficient and denoted by k. This 
is a quantity of phenomenological interest. Measurements at the Relativistic Heavy Ion Collider 
(RHIC) have shown that heavy quarks display substantial elliptic flow [l[ 3], implying stronger 
medium interactions than extrapolated weak-coupling calculations would suggest. It has been 
estimated [3] that a diffusion coefficient D < 1/T, or equivalently, k/T 3 > 2, is required in order 
to accomodate the RHIC data. Theoretically, we will be concerned with collisional energy loss, 
which is expected to dominate for small P/M; however one should keep in mind that in heavy- ion 
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collisions, the quarks are not produced at rest, hence radiative energy loss can play a role too 



see [4] for a recent review. A separate motivation is that the dynamics of the heavy-quark probe is 
somewhat simpler to study than the diffusion of transverse momentum carried by the constituents 
of the plasma, which is characterized by the shear viscosity. There are reasons to expect that 
once the heavy quark diffusion constant has been determined, it can be used to calibrate other 
transport coefficients such as the shear viscosity, because ratios of transport coefficients are more 
stable predictions of the weak-coupling expansion . 

We briefly review the known analytic results for the momentum diffusion coefficient. To strict 
leading order, k is given by [fj (Cf = %\?—) 

g 2 C F T 2 / 2T 1 C'(2) a iVflog_2 \ m 
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Unfortunately, for realistic values of the Debye screening mass mo, this expression is negative. 
Carrying out the integrals without expanding inmo { m % = 9 2 T 2 (N c /3 + iVf/6) in leading order), 
positive values are obtained, for instance at g 2 ~ 3, k/T 3 ~ 0.4 for JVf = 3 and 0.23 for Nf = 5|. 
Remarkably, a next-to- leading order computation has been carried out [6|, with the result k/T 3 ~ 
2.7 for g 2 ~ 3.0 for Nf = 3. It is unclear at this point how useful the expansion is. 

The N = 4 SYM theory provides an interesting testing ground for analytic methods. The 
result obtained by holographic methods in the strong coupling limit of N = 4 SYM theory is 
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a], where A = g 2 N c is the 't Hooft coupling. On the other hand, the NLO 



weak-coupling result has also been worked out 



A 2 T 3 / 1 \ 
SYM = log + 0.4304 + 0.8010VA 
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(2) 



Although the apparent convergence is again poor, it is encouraging that the strong coupling result 
crosses the weak-coupling result at an intermediate coupling of A ~ 10.2. This suggests that the 
values obtained from NLO expressions such as ([2]) are in the right ball-parJ3. 

In the next section, we review the steps that lead to a Kubo formula for the momentum diffusion 
coefficient k in the static limit of the probe-quark. In section III we present our numerical results 
for the Euclidean correlators obtained by Monte-Carlo simulations, and conclude in section IV. 



1 A ~ 10.2 is a moderately large coupling. For a 't Hooft coupling of 10-12, the spatial scalar correlator computed 
by holographic methods in the SYM theory and on the lattice in the SU(iV c ) gauge theory exhibits a non-trivial 
agreement up to about 2T C [101 ]. In this range, the 0(1/A 3 ^ 2 ) correction to the A = 00 shear viscosity to entropy 
density ratio amounts to about 20%. 
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II. THE SPECTRAL FUNCTION OF THE HEAVY-QUARK CURRENT 

This section is mainly a review of recent literature on the subject of heavy-quark diffusion, with 
some comments pertinent to calculations performed in Euclidean space. The basic definitions of 
thermal correlators and the relations among them are gathered in the Appendix, to which we refer 
the reader for unexplained notation. The main equations of the linear response framework are also 
reviewed there. 

To study the diffusion of a conserved charge such as the heavy quark number N = J dx n(t, x), 
the Hamiltonian can be perturbed by 

H tl = H-J 6x/i(t,x)n(t,x), fj,(t, x) = n(x) e et 6(-t). (3) 

In the hydrodynamic treatment of the problem, Fick's law j = -DVn and the conservation 
equation dtn + V • n = lead to the diffusion equation, whose solution in Fourier space takes the 
form 



Xs(k)n(k) |„\ — f°° Jut f J_„-tfc-iC 

-iui + Dk 2 



n(u,k)= , , „/ 2 , n(oj,k)= j e iujt j dxe- ik - x n(t,x). (4) 



o 

Here Xs(k) = ft f dx e~ lkx (n(t,x)n(0)} and Xs = Xs(0) is the particle number susceptibility. Via 
(|A20p . this determines the retarded correlator for small oj and k, 

„ nnf (Dk 2 ) 2 + iu Dk 2 

Before proceeding further, it is also instructive to write down the correlator in the real-time domain, 

G^(t,k) 1 ~°° X s{k)Dk 2 exp ( - DkH). (6) 



Using Eq. (|A10j) . the contribution of this exponential tail to the Euclidean correlator (say) at 
t = ft/2 is Xs(k), independently of the diffusion coefficient D. This is a manifestation of the 
difficulty to extract transport information from Euclidean correlators. 

The longitudinal part of the current correlator (i.e.{j z j z ) if k = (0, 0, k)) is related to the density 
correlator by the current conservation equation, pl(oj, k) = ^p nn (u, k) = ImG^ n (cj, k). Thus 
the current spectral function reads 

p L (u,k) Xs{k) Dlj 2 



W 7T Ujl + (Dk 2 ) 2 ' 

implying in particular the Kubo formula 



(7) 



D Xs = vr lim lim pL ^^ k \ ( 8 ) 
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The diffusion of a heavy quark (one has in mind the charm or preferably, from a theoretical 
point of view, the bottom quark) in the quark-gluon plasma is characterized by a time scale M/T 2 
which is long compared to the thermal time scale of 1/T. For this reason it is expected that a 
classical Langevin equation should appropriately describe the thermalization of heavy quarks Qj. 
See and References therein for a derivation. The heavy-quark's classical equations of motion 



arc 



dx p dp J _. . . . 

(«W(*')> = -f)- (io) 

For a given £ (i) , the equation is easily solved to give 

p{t) = e-* t [p(Q) + jl&si{s)e* s ], (11) 

implying 

lim{ Pi (t)p J (t)) = ^5 ij . (12) 

The equipartition of energy requires to be |T in equilibrium; the drag and fluctuation coeffi- 
cients are thus related by the fluctuation-dissipation relation 

(13) 
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The mean square distance covered by the particle is also easily worked out. For a thermal initial 
distribution of momenta, (pi(0)pj(0)) = MT5ij, it reads 

\{x\t)) = 2D[i-±(l- e"*)]. (14) 

This equation describes both the early-time directed motion, ^(x 2 (t)) = \v 2 t 2 , ^v 2 = jj, and the 
late-time diffusive motion, ^(x 2 (t)) = 2Dt [3]. Let now P(t,x) be the probability that a heavy 
quark starts at the origin at t = and moves a distance x over a time t. If the distribution of 
heavy quarks at time zero is n(0, x), at time t it will be given by the convolution 

n(t,x) = J dx'P(t,x - x')n(0,x'), (15) 

or equivalently 



n(t,k) = P(t,k)n(0,k). 



(16) 
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If one assumes the noise to be Gaussian distributed, then the probability distribution P(t, x) is 
Gaussian [l^ . with a width given by Eq. (|14p . and therefore so is P(t,k). Applying the general 
rule (|A20p . one then obtains the retarded correlator, 



1 + icj dt e lujt P(t, k) 
Jo 



(17) 



The resulting spectral functions were obtained numerically in [12]. In particular, the k = spectral 
function takes the form 



p L (u,0) _ Xs T rj 



(18) 
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The spectral structure that is obtained from the Langevin equation is expected to arise for a 
sufficiently heavy diffusing particle. Conversely, the presence of a transport peak allows one to 
define the quantities appearing in the Langevin equation directly from the spectral function [l^ ]. 
The effective mean-square velocity is then given by 

^ 2 > - f r ^plh, (19) 

Xs J -A w 

where A is a cutoff that separates the scale rj from the correlation time of the medium (which is 
typically of order T, or gT at weak coupling). The 'kinetic mass' M^ n is further defined so as 
to satisfy the equipartition theorem, M^[ n (v 2 ) = 3T. Finally, the momentum diffusion coefficient 
k{M) can be defined as 

2ttM?- 

k{M) = ^ cop L (co) . (20) 

Xs r?<M<A 

A weak-coupling calculation shows that while k(M) and D are only weakly dependent on M, the 
drag coefficient rj ~ T 2 /Mx a power of the coupling constant is parametrically small compared to 
the medium time-scale. By reexpressing the right-hand side of (1201 ) in terms of the correlator of 

nn 

two heavy-quark currents, the authors of [31,113] were able to formulate the task of computing n in 
the static limit of Heavy-Quark Effective Theory (HQET). In that limit, the leading contribution 
is the 'force-force' correlator, where the force is given by the Lorentz expression gE. The relevant 
Euclidean correlator reads, after evaluating the fermion line contractions, 



( |( )rr x Re Tr {U(J3, t)gE k (t, 0)U(t, 0)gE k (0, 0)) / 
° E {t) ~ -3 <ReTrtf(/3,0)) ' (21) 

where the color parallel transporters U(t2,t\) in the fundamental representation are propagators of 
static quarks. In particular the Polyakov loop appears in the denominator of (|2ip . The momentum 
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diffusion coefficient is given by the low-frequency limit of the corresponding spectral function via 
Eq. flSHD , 

K = li m ^ p HQET (w) _ (22) 

The obvious advantage of this formulation is that the large scale M has disappeared from the 
problem. The spectral function p H Q ET has been studied in detail at next-to-leading order (NLO) 
in perturbation theory [l^j]. Remarkably, even in the weak-coupling limit, the function is smooth 
as small frequencies. This is in contrast with the narrow transport peaks that are found at weak- 
coupling in e.g. the shear channel. This property represents a clear advantage for numerical studies 
of the spectral function: the form of the kernel in Eq. (|A12|) makes the Euclidean correlator very 
insensitive to the functional form of the spectral function at u < T. On the downside, while 
the spectral function of the current-current correlator grows as j^z^ 2 at large frequencies, the 
spectral function of the S-field correlator grows even faster, like ^§fuJ 3 Q. This implies that the 
low-frequency part makes a comparatively small contribution to the Euclidean correlator studied 
in the next section. 



III. FORCE-FORCE CORRELATORS FROM LATTICE GAUGE THEORY 



In this section we describe a first calculation of the HQET 'force-force' correlator (|2ip by means 
of Monte-Carlo simulations in Euclidean space. The calculation is performed in the deconfined 
phase of SU(3) gauge theory. We employ the isotropic Wilson action [151 ]. 



S g = iE Tr{l-P^(x)} 



.95 



(23) 



where the 'plaquette' P^ u is the product of four link variables U^{x) £ SU(3) around an elementary 
cell in the (^,^) plane. The size of the lattice is N T x N%, with periodic boundary conditions 
in all directions. As a local update algorithm, we use the standard combination of heatbath and 
over-relaxation 1614191] sweeps in a ratio increasing from 5 to 7 as the lattice spacing is decreased. 



No multi-level algorithm 



20, 



211 ] was used here, although we expect that for sufficiently small 



lattice spacing, such an algorithm will be beneficial, since the fluctuations of the E fields will be 
UV-dominated. To set the scale we use the parametrization of the Sommer scale ro ~ 0.5 fm g iven 



m 



221 ] based on the data [23J, and convert between ro and T c by using vqT c = 0.746(7) ([24| and 



Refs. therein). 
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FIG. 1. The Euclidean force- force correlator extracted from simulations with the Wilson action at /3 = 7.483 
and the HYP1 discretization of HQET. The black curve indicates the NLO result of Burnier et al [14]. This 
result has been used to determine the normalization of the E-field by matching the lattice data to it at 
tT = j and 3.1T C . Statistical error bars are smaller than the data symbols. 



There is a lot of freedom in discretizing the HQET Lagrangian (see 



2511 for an introduction). As 



color parallel transporters, we choose the lattice gauge links obtained after one iteration of HYP 



smearing [261 ] . This action has been used extensively at zero-temperature [27|, and our choice of 
smearing parameters is known in the literature as the HYP1 action [2a |. The original motivation 
for choosing this action is that it reduces the size of the UV-divergent self-energy, thus leading 
to a reduction in the statistical errors. At the same time the cutoff effects were found to be 
controllable. In 



291 ] . the chromo-magnetic field B was also discretized with HYP-link variables, 
yielding also a benefit in statistical error reduction. Following this example, we use HYP-links for 
the chromo-electric field as well. 



While the renormalization factor for the B field was computed in [29|] (it is scale-dependent), the 
corresponding factor for the chromo-electric field remains to be computed. The same methods apply 
to the renormalization of this factor, but in the mean time we use a preliminary way of normalizing 
the chromo-electric operator. We expect that at short time separations t, perturbation theory 
provides an accurate prediction for Ge(i). Therefore the absolute normalization of the chromo- 
electric field can be obtained by requiring that the lattice correlator match the NLO perturbative 
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FIG. 2. The Euclidean force-force correlator, normalized as in Fig. [TJ at t — 1/2T, with /3 = 7.483 and 
N T = 8, 12, 16, 22. Statistical error bars are smaller than the data symbols. 



prediction at some reference time t ie {. The difficulty in this procedure is that i re f must still be 
large enough in lattice units for the discretization errors to be under control. In this work we made 
the compromise to choose i re f = 1/4T at a temperature of about 3.1T C on a 16 x 64 3 lattice. In 
physical units, this represents a separation of about 0.05fm, and four lattice spacings in lattice 
units. At one value of the bare coupling g^, this normalization is then valid for other values of iV T . 
In this way, we have varied the temperature by changing N T from 8 to 22 at fixed 6/<?q = 7.483. 
The chromo-electric field correlators are displayed in Fig. [TJ Only those points with t/a > 4 
are displayed. At shorter separations, the lattice correlators exhibit non-monotonicity (and are 
even negative at separation and 1 lattice spacing), a not unexpected lattice artefact. It is clear 
that with the procedure adopted here, we cannot disentangle discretization errors from the bare- 
coupling dependence of the renormalization factor of E. Our results should accordingly be regarded 
as preliminary. Eventually the renormalization factor should be computed along the lines of 
where the B operator was treated instead. 



29], 



To repeat, the data at 3.1T C has been calibrated to match the 0(g 4 ) result of [141 ] at t = 1/4T. 
In Fig. [JJ different temperatures are accessed by varying N T . On the logarithmic scale of the plot, 
the temperature dependence of the Euclidean correlator normalized by T 4 is weak. Furthermore 
the NLO perturbative prediction provides a rather good description of the i-dependence at 3.1T C . 
The temperature dependence of the mid-lattice data points (normalized as in Fig. [JJ) is shown in 
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FIG. 3. The quantity displayed is Q(t), denned in Eq. The '+' denote the 0(.g 4 ) prediction at 

3T C [l^j]. The 'x' denote the function fl(t) that results from adding the low- frequency contribution Ap(cj) = 
i-Are tanh(^)#(A — \u\) to the 0(<7 4 ) spectral function, with parameters An and A given in the caption. 
Statistical error bars are smaller than the data symbols, but discretization errors are probably non-negligible, 
particularly at the smaller t values. 

much greater detail in Fig. [2j The temperature variation between 2T C and 6T C is on the order 
of 20%. The sign of the variation is the same as predicted by the perturbative expression: the 
magnitude of the force-force correlator in the gluon plasma increases as the temperature decreases. 

Although we do not yet have non-perturbative control over the absolute normalization of the 
correlators, the latter cancels out in the relative fall-off of the correlator. An observable that 
measures this is the quantity £l(t) > defined by 

G E {t - a/2) _ cosh [«(t)09/2 - (t - a/2))] 



G E (t + a/2) cosh [Q(t)(/3/2 - (t + a/2))] 



(24) 



We remark that Q(t) has a continuum limit, in which Q tanh S~2(/3/2 — t) = — ^ log G^;(t). One can 
interpret it as the location of a delta function in the spectral function which by itself reproduces the 
local fall-off of the Euclidean correlator. The function f2(i) is displayed in Fig. [3] for three different 
temperatures, where the temperature is varied this time by changing the bare coupling g$ at fixed 
iV T = 16. We note that because the statistical samples of the numerator and the denominator 
in Eq. (|24p are highly correlated, the numerical results for these ratios have uncertainties at the 
few-permille level. Figure [3] also displays the perturbative prediction at 3.1T C corresponding to the 
curve appearing in Fig. [TJ it is based directly on Eq. (|24p rather than on the continuum version 
of this equation, to allow for a more direct comparison with the lattice data. The lattice £l(t) 
differs from the perturbative one only by a few percent, namely the latter falls off slightly more 
steeply. Since the difference is so small, it could partly be due to discretization effects. To reduce 
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their influence, we concentrate on the largest values of t. One may ask nonetheless, how large a 
difference in the transport coefficient k could this discrepancy possibly correspond to. An estimate 
is obtained by adding a low-frequency correction to the perturbative spectral function, 

Ap(u) = i • Ak tanh(^) 0(A- |w|). (25) 

Other functional forms than (|25p . such as a Breit-Wigner curve, would perhaps be more realistic, 
but would not change our conclusions in any significant way. Adding such a term to the spectral 
function has the effect of making the Euclidean correlator flatter, and indeed, by adjusting Ak, one 
can obtain good agreement for the largest two t values between the perturbative prediction modified 
by Eq. fl25J) and the lattice data. At T = 3.1T C we find that, for A = T and Ak/T 3 = 0.352(38), 
agreement is obtained with the lattice data at the two largest t values. This represents a substantial 
enhancement of k over the leading-order perturbative value mentioned in the introduction. An 
equally good agreement is obtained if one chooses A = 2T, which leads to Ak/T 3 rj 0.204(22). 



IV. CONCLUSION 



To summarize, we have found that the Euclidean force- force correlator evaluated at t = (3/2 
admits a « 20% increase as the temperature is lowered from 6 to 2T C . At 3.1T C its t-dependence is 
described at the few-percent level by the recent NLO perturbative result [14J]. The somewhat flatter 
behavior seen in the lattice data can be explained by an enhancement of the spectral function at 
low frequencies, and adopting this explanation leads to a substantial increase of k over the leading- 
order result obtained in ^]. While it is too early to draw phenomenological conclusions, at present 
the increase appears to be not quite sufficient to explain the experimentally observed elliptic flow of 
heavy quarks, as discussed in the introduction. It will be interesting to see the results brought by 



30] 



the current lead-lead collisions at the LHC. The ALICE experiment has very recently reported 
an integrated elliptic flow (i.e. of light constituents) about 30% larger than at RHIC. 

From a technical point of view, it is encouraging that a precision at the few per-mille level can 
be achieved on the force-force correlator on a 16 x 64 3 lattice. Given this statistical precision, 
it is essential, as a next step, to control the discretization errors at a comparable level. One 
could then study more systematically and quantitatively the implications of the lattice data for 
the spectral function. A mandatory step in controlling the discretization errors is to compute 
the renormalization factor of the E field non-perturbatively. The motivation to carry out these 
calculations is now quite strong. 
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Appendix A: Thermal correlators and linear response 

We start by recalling some of the definitions, which allows us to fix our notation. At finite 
temperature T = 1/(3, correlation functions are defined as 

G AB (t) = Tr{pA(t)B(0)}, (Al) 

with p = the equilibrium density matrix. Expectation values of commutators, 

G AB (t) = m{p[A(t),B(0)}} = * (G> B (t) ~ G BA (-t)) , (A2) 

play a particularly important role in finite-temperature physics. The integral transform over the 
positive half-axis 

POO 

G A P{u)= / dte iujt G AB (t) (A3) 
Jo 

is analytic in the half complex plane Im(w) > 0. We will refer to it as the retarded correlator. The 
spectral function, which is really a distribution, is defined as 



1 



+oo 



p ab (uj) = / dt e iujt G AB (t). (A4) 

For B = , the spectral function is identically related to the imaginary part of the retarded 
correlator, 

P AA1 ( w ) = I Im Gjf (w)el. (A5) 
The Euclidean correlator is defined as 

G AB (t) = G AB (-it). (A6) 
It obeys the Kubo-Martin-Schwinger relation 

G BA ((3-t) = G AB (t), (A7) 
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and therefore admits a representation as a Fourier series, 

G E {t) =tJ2 g e e ' iuJet > g e ] = / dt e luJet G E (t) , (A8) 

where W£ = 2nT £ and we have dropped the label specifying the operators A, B. 
In frequency space, the Euclidean and retarded correlators are related by 

G R (kJt) = Gf • (A9) 

The I = case has to be treated somewhat more carefully, see [3l| for details. In coordinate space, 
the relation between the Euclidean correlator and the real-time correlator is 

an*) + os*w - rjy n + G *v,) CMh(2 ^, ( ) 2 ::2 (2jrt) . (mo) 

Haif {t) - crw) - r jT d W) .^ff^W <ah) 



Finally, the most commonly used relation between Euclidean and real-time correlators is the mixed 
coordinate-frequency space relation 

<*»<«) + - «) » ^|t^. (a*> 

G - (t) _ G * BW _ t) =£y » fi^|_!> . (A13) 

The retarded correlator G^ B is important because it is related to the response of operator A to 
a perturbation of the system by operator B. A time-dependent perturbation of the Hamiltonian 
by an operator B, 

H f (t) = H - f(t)B(t), (A14) 

leads to a 'response' of physical quantities, i.e. a change in their expectation values with respect 
to the unperturbed ensemble. The evolution equation of an operator A is given by 

i^A(t) = -[H f (t),A(t)]. (A15) 

One then finds that to linear order in /, the expectation value of A in the perturbed system is 

6(A(t)) = (A(t)) f - (A(0)> = f dt'G AB (t - t')f(t') + 0(f). (A16) 

J — oo 

Equation (|A16|) is the master formula of linear response theory. It shows that the retarded cor- 
relator Gft B determines the response of an observable A to a time-dependent external field that 
couples to B. A source term of the form 

f(t) = e^6(-t)f (A17) 
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is often adopted to study how the system relaxes back to equilibrium after having been perturbed 
adiabatically. The static susceptibility is defined as the expectation value of A at t = 0, 

S(A(t = 0)) f = xt B fo. (A18) 

From (|Al"6j) . it follows that 

/•oo 

X f B = / dte-^G AB (t) = G% B (ie). (A19) 
Jo 

Integrating both sides of (|A16p . / °° dwe iIj( (.), one obtains for the adiabatic perturbation ()A17|) 

poo 

Gi B (uj)f = {6A(0))f +iu dte iujt (5A(t)) f . (A20) 



o 

This formula shows that the relaxation of observable A back to its equilibrium value and the 
retarded correlator G^ B (u) are in one-to-one correspondence. Since the late-time relaxation is 
expected to be described by hydrodynamic evolution, this equation can be exploited to obtain a 
prediction of the small-w functional form of Gr. 
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